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fSj Abstract. In the current paper, we further improved the model for the migration of 
planets introduced in Del Popolo et al. (2001) and extended to time-dependent planetes- 
■ imal accretion disks in Del Popolo and Ek§i (2002). In the current study, the assumption 
of Del Popolo and Ek§i (2002), that the surface density in planetesimals is proportional 
'" p^ ^1 to that of gas, is released. In order to obtain the evolution of planetesimal density, we 
Q use a method developed in Stepinski and Valageas (1997) which is able to simultaneously 

' coupled to migration model introduced in Del Popolo et al. (2001) in order to obtain the 



follow the evolution of gas and solid particles for up to lO^yrs. Then, the disk model is 



X 



migration rate of the planet in the planetesimal. We find that the properties of solids 
known to exist in protoplanetary systems, together with reasonable density profiles for 
the disk, lead to a characteristic radius in the range 0.03 — 0.2 AU for the final semi- major 
axis of the giant planet. 
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1. Introduction 

The discovery of solar-like stars showing evidences for planets orbiting around 
them dMayor fc Queloz 1995||Marcy et al. 2000[|T^et al. 2000||Butler et al. 2001|l has 

greatly intensified the interest in understanding the formation and evolution of planetary 
systems, as well as the long-standing problem of the solar system origin. 
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The extra-solar planets discovered so far are all more massive than Saturn, and most 
either orbit very close to their stars or travel on much more eccentric paths than any of 
the major planets in our Solar System. 

It is difhcult to explain the properties of these planets using the standard 
model for planet formation l|Lissauer 19981 IBoss 1995jl . Current theories (jMizuno 19801 
IBodenheimer fc Pollack 1986|l predict that giant planets were formed by gas accretion 
onto massive (~ 15Mq) rocky cores which themselves were the result of the accumulation 
of a large number of icy planetesimals. The most favourable conditions for this process 
are found beyond the so-called "snow line" ( |Hayashi 1981| ISasselov fc Lecar 2000|l . As 
a consequence, this standard model predicts nearly circular planetary orbits and giant 
planets distances > 1 AU from the central star where the temperature in the protostellar 
nebula is low enough for icy materials to condense (|Boss 19951 IBoss 19961 but see also 
IWiichterl 1998TlWiichterl 199fill. 

Therefore, in the case of close-in giants, it is very unlikely that such planets were 
formed at their present locations. Then, the most natural explanation for this paradox, 
and for planets on very short orbits, is that these planets have formed further away in the 
protoplanetary nebula and they have migrated afterwards to the small orbital distances 
at which they are observed (see DPI and DP2 for a detailed discussion of migration 
mechanisms). In particular, in DPI and DP2 we showed that dynamical friction between 
a planet and a planetesimals disk is an important mechanism for planet migration and 
we pointed out that some advantages of the model are: 

a) Planet halt is naturally provided by the model. 

b) It can explain planets found at heliocentric distances of > 0.03 — 0.04 AU, or planets 
having larger values of eccentricity. 

c) It can explain metallicity enhancements observed in stars having planets in short- 
period orbits. 

d) Radial migration is possible with moderate masses of planetesimal disks, in contrast 
with other models. 



In DPI, following Opik 1976 it was assumed that the surface density in planetesimals 
Ss varies as Ss(?') — 5]o(lAU/r)'^/^, where E©, is the surface density at 1 AU. In DP2 
the previous assumption was substituted by a more reliable model for the disk, and in 
particular we used a time-dependent accretion disk, since it is widely accepted that the 
solar system at early phases in its evolution is well described by this kind of structure. 
An important assumption of DP2 was that the surface density in planetesimals remains 
proportional to that of gas: T^s{R, t) oc t). However, it is well-known that the distri- 
bution of planetesimals emerging from a turbulent disk does not necessarily reflect that 
of gas (e.g., |Stepinski fc Valageas 1996[ |Stepinski fc Valageas 1997| ). Indeed, in addition 
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to gas-solid coupling, the evolution of the distribution of solids is also godetermined by 
coagulation, sedimentation and evaporation/condensation. In order to take into account 



to simultaneously follow the evolution of gas and solid particles for up to lO^yr. The 
main approximation used in this model is to associate one grain size to a given radius 
and time. Then, we use the radial distribution of planetesimals given by this model to 
evaluate the planet migration, which is calculated as in |Del Popolo et al. 2001| 

This paper is organized as follows. In Sect. 2, we describe the disk model we use 
to obtain the distribution of the planetesimal. Then, in Sect. 3, we briefly review the 
migration model introduced in [Del Popolo et al. 200 1| Finally, we describe our results in 
Sect. 4 and Sect. 5 is devoted to conclusions. 

2. Disk model 

It is well-known that protostellar disks around young stellar objects are common: between 
25% to 75% of young stellar objects in the Orion nebula seem to have disks with mass 
IQ-^M© < Afd < lO^^M© and size 40 ± 20 AU |Beckwith fc Sargent 1996| Moreover, 
observations of circumstellar disks surrounding T Tauri stars support the view of disks 
having a limited life-span and characterized by continuous changes during their life. These 
evidences have led to a large consensus about the nebular origin of the Solar System. 
Moreover, it clearly appears necessary to model both the spatial and temporal changes 
of the disk (which cannot be handled by the minimum-mass model nor by steady-state 
models). Besides, one also needs to describe the global evolution of the solid material 
which constitutes, together with the gas, the protoplanetary disk. ^ As usual, the time 
evolution of the surface density of the gas S is given by the familiar equation (e.g., 
IStepinski fc Valageas 1997| ): 



where ft is the turbulent viscosity. Since vt is not an explicit function of time, but instead 
depends only on the local disk quantities, it can be expressed as vt = ft(S, r) and Eq.(^ 
can be solved subject to boundary conditions on the inner and outer edges of the disk. 
The opacity law needed to compute i^t is obtained from IRuden fc Pollack 19911 Then, 
Eq.QJ is solved by means of an implicit scheme. Note that the evolution of the gas is 
computed independently from the evolution of particles (which only make 1% of the 
gas mass). Next, from S(r, i) we can algebraically find all other gas disk variables. 

^ The knowledge of this distribution and its time evolution is important to understand how 
planets form and in this paper it is a key issue since we wish to study the planet migration due 
to the interaction between planets and the local distribution of solid matter. 





(1) 
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Next, as described in |Stepinski fc Valageas 1997| the evolution of the surface density 
of solid particles Eg is given by: 



dt r dr 



or 



1 d 




r dr 





(2) 



The first diffusive term is similar to Eq. JQ) , where the effective viscosity z/g is given by: 



/ 9 



= with sc = (f + r!kts) y 1 + (3) 

Here we introduced the Schmidt number Sc which measures the coupling of the dust to 
the gas turbulence. We also used the relative velocity v between a particle and the gas, 
the turbulent velocity Vt, the Keplerian angular velocity ilk and the so-called stopping 
time ts- The average (..)s refers to the vertical averaging over the disk height weighted 
by the solid density. The dimensionless quantity (flk^s) measures the coupling of the 
solid particles to the gas. Therefore, the evolution of the dust radial distribution can be 
significantly different from the behaviour of the gas, depending on the particle size (see 
IStepinski fc Valageas 1996| for a detailed study). 

The second advective term in Eq. Q is due to the lack of pressure support for the dust 
disk as compared with the gas disk. Thus, it is proportional to the azimuthal velocity dif- 
ference U0 between the dust and the gas. We refer the reader to |Stepinski fc Valageas 19971 
for a more detailed presentation, see also lKornet et al. 20011 

In this article we are mainly interested in the distribution of solids at small radii 
hence we consider only one species of solid particles: high-temperature silicates with 
Tcvap = 1350 K and a bulk density Pbuik = 3.3 g cm~'^. Thus, in our simplified model 
we follow the evolution of three distinct fluids: the gas, the vapour of silicates and the 
solid particles. In this way, we obtain the radial distribution of the planetesimal swarm 
after 10^ yr. Of course, at these late times when planetesimals have typically reached 
a size of a few km or larger, gravitational interactions should play a dominant role 
with respect to coagulation. However, if these interactions do not significantly affect the 
radial distribution of solids (note that the radial velocity of such large particles due to the 
interaction with the gas is negligible) we can still use the outcome of the fluid approach 
described above to study the migration of giant planets, as detailed below. 



3. Migration model 

In order to study the migration of giant planets, the model developed in DPI is used 
(see also DP2). Since this model has already been described in these two papers, we 
only recall here the main points. We consider a planet revolving around a star of mass 
Af, = IMq. The equation of motion of the planet can be written as: 

r = Fo + R (4) 
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where the term Fq represents the force per unit mass from the star, while R is the 
dissipative force (i.e. the dynamical friction term-see IMelita fc Woolfson 1996)1 . 

In order to take into account dynamical friction in a disk structure, we use Binney 
(1977) model: 

K ^ -k\ivnie\i ~ k±vi±e± (5) 

where ey and ej^ are two unit vectors parallel and perpendicular to the disk plane and 
fc|| and k± are given in Binney (1977) and DPI. 

Since the damping of eccentricity and inclination is more rapid than radial migration 
Pda 19901 Hda fc Makino 19921 DPI), we only deal with radial migration and we assume 
that the planet has negligible inclination and eccentricity {ip Cp ~ 0) and that the 
initial distance to the star of the planet is 5.2 AU. We do not need to follow the evolution 
of the size distribution of planetesimals since for m <^ M the frictional drag obtained 
does not depend on the mass m of the planetesimals, because the velocity dispersion only 
depends on the mass M of the giant planet. We merely use the planetesimal density ps 
reached after 10^ yrs, assuming that the height of the planetesimal disk does not evolve 
significantly. 

4. Results 

In this article, similarly to DPI and DP2, we are mainly interested in studying the 
planet migration due to the interaction with planetesimals. Our model starts with a fully 
formed gaseuos giant planet of 1 Mj at 5.2 AU. For this reason we assume that the gas 
is almost dissipated when the planet starts its migration. ^ While the gas tends to be 
dissipated, (several evidences show that the disk lifetimes range from 10^ yr to 10^ yr, see 
IStrom et al. 1993linuHen fc Pollack 1991|l . the coagulation process induces an increase of 
the density of solid particles with time (see Fig.^ and gives rise to objects of increasing 
dimensions. 

Once solids are in the form of planetesimals, the gas coupling becomes unimportant 

and the radial distribution of solids does not change any more. This is why we do not need 

to calculate its evolution for times longer than 10^ years. Then the disk is populated 

by residual planetesimals for a longer period. Here it is important to emphasize that 

planetesimal formation is not independent from initial conditions. In particular, the final 

solid surface density depends in an intricate way on the initial disk mass Md and on the 

turbulent viscosity parameter a, (see |Stepinski fc Valageas 1997| and lKornet et al. 200l)l . 

^ Clearly the effect of the presence of gas should be that of accelerating the loss of angular 

momentum of the planet and to reduce the migration time. 
^ Note however that the size distribution of planetesimals keeps changing due to their mutual 

gravitational interaction. 
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Table 1. Properties of the initial gas disk 
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In order to investigate the dependence of the giant planet migration on the properties 
of the protoplanetary disk we integrated the model introduced in the previous sections 
for several values of the initial disk surface density (i.e. several disk masses), and different 
values of a. More precisely, as in [Stepinski fc Valageas 1997| we consider an initial gas 
surface density of the form: 

So(r) = Si [1 + (r/ri)2] + ^^(r/lAU)-!-^ (6) 

The quantities Si, ri and S2 are free parameters which we vary in order to study different 
disk masses. The values used are given in Tab. 1, where Mfi is the gas disk mass (in units 
of Mq), J50 is the disk angular momentum (in units of 10^° g cm^s-^), Si and S2 are in g 
cm-^ and ri is in AU. The first term in Eq.® ensures that there is some mass up to large 
distances from the star, while the second term corresponds to the central concentration 
of the mass and determine the location of the evaporation radius. Note however that in 
any case the evaporation radius for the high-temperature silicates we study here remains 
of order 0.1 AU. As previoulsy explained, we consider only one species of solids: high- 
temperature silicates with T^vap = 1350 K. We initialize the dust subdisk at time t ~ 10^ 
yrs (i.e. after the gas distribution has relaxed towards a quasi-stationary state) by setting 
the solid surface density Sg as: Sg = 6 x 10-^S in order to account for cosmic abundance. 

We show in Fig. ^ the evolution of the solid midplane density for the case Md = 
1O-^M0, and for different values of a. We can see a converged radial distribution of 
solids emerge at late times of order 10^ yr. Note that although the radial distribution 
of planetesimals depends on the value of a, the total mass of solids in a disk is roughly 
independent of a since it remains approximately equal to the initial mass of solids in 
the disk. This means that solids are reshuffled within the disk but they are not lost 
into the star. ^ The value of a determines the radial distribution of solids: particles in a 
disk with a larger value of a (more turbulent disk) have larger inward radial velocities 
and consequently are locked into planetesimals closer to the star than particles in a less 
turbulent disk. Thus, the smaller the value of a the broader the final distribution of solids. 

* In fact, solids initially located close to the evaporation radius are lost but they constitute a 
small percentage of the total solid material, which is predominantly located in the outer disk. 
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Fig. 1. (a) Evolution of solids for a disk with Md = O.IM© and a = 0.1 a.t t = lO'^yrs 
(dashed line), t = lO^yrs (dotted line) and t = lO^yrs (solid line), (b) Same as Fig. la 
but with a — 0.01. (c) Same as Fig. la but with a ~ 0.001. (d) Same as Fig. la but with 
a = 0.0001 



The evaporation radius is located at ~ 0.1 AU while the outer limit of converged Ps is 
about 10 AU (it moves outward for smaller values of a, (e.g., ~ 70 AU for a = 10^*)). 

The coagulation process gives rise to solids of 10^ — 10^ cm. In disks characterized by 
smaller values of a, and thus a more extended distribution of solids, the range of sizes goes 
from 10^ — 10^ cm at the evaporation radius down to lO'^ — 10^ cm at the outer limit. This 
is because the coagulation process is less efficient at larger radii where the solid density is 
smaller and the velocity dispersion of the dust decreases (along with the gas temperature 
which governs the turbulent velocity). At later times these solids will continue to increase 
their sizes, but will not change their radial position, as they are already large enough to 
have a negligible radial motion. We refer the reader to [Stepinski fc Valageas 1997| and 
IKornet et al. 2001l for more detailed discussions of the behaviour of protoplanetary disks. 

Now, using the converged radial distribution of planetesimals derived in the previous 
section, we show in Fig. |21the evolution of the semi-major axis a{t) of a 1 Mj planet 
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in a disk for a = 0.1 and Aid = 0.1 (solid line), 0.01 (dotted line), 0.001 (short-dashed 
line), O.OOOIM0 (long-dashed line), respectively. We recall here that the planet is initially 
located at 5.2 AU with ip ~ ep ^ 0. 

It can be clearly seen that for a fixed value of a a disk of larger mass leads to a more 
rapid migration of the planet. This behaviour is quite natural since a more massive disk 
obviously yields a stronger frictional drag. Thus, in the cases M = 0.1, O.OIA/q the planet 
migrates to ~ 0.08 AU in ~ 1.5 x 10^ yr and to ~ 0.03 AU in ~ 2.5 x 10^ yr, respec- 
tively. When the planet arrives at this distance the dynamical friction switches off and 
its migration stops. The stopping is simply due to the inner radius of the planetesimal 
disk. The latter is set by the evaporation radius. Indeed, solid bodies cannot condense 
at such small orbital radii r ^ 0.1 AU because the temperature is too high. Of course, 
this evaporation radius Tcvap depends on the properties of the solid grains we comsider. 
For instance, for ice particules we have rovap ~ 1 AU (e.g., |Stepinski fc Valageas 1997| 
IKornet et al. 200T)l . In this article we wish to understand the small orbital radii of ob- 
served planets, over the range 0.03 — 0.15 AU. Therefore, we are interested in the inner 
regions of the disk where only high-temperature silicates with Tcvap ^ 1350 K survive. 
This is why we selected this component in this study. Then, the main point of Fig.|21is 
that the properties of solids known to exist in protoplanetary systems, together with rea- 
sonable density profiles for the disk, lead to a characteristic radius in the range 0.03 — 0.2 
AU for the final semi-major axis of the giant planet. Note in particular that this process 
naturally explains why the migration stops at such radii. 

For less massive disks, Md = 10^^, IO^^Mq, the migration is slower and the planet 
has not enough time to migrate belows ~ 4 AU. Note that we stop our calculation 
after 4 x 10^ yr which is the typical age of protoplanetary disks like ours (the Sun is 
^ 4.5 X 10^ yr old). Moreover, the planetesimal disk should have cleared out by this 
time. It is interesting to note that the planet moves closer to the central star in the case 
Al = Q.OIAIq than in the case Al = O.lAf0. This is due to the fact that less massive disks 
usually have a smaller surface density which leads to a smaller temperature. This in turn 
implies a smaller evaporation radius so that the planet can move closer to the star. ^ In 
order to study the effect of viscosity on migration, we performed three other calculations 
with a = 10^^, 10^'^ and 10^*, also plotted in Fig. |21 We can note that the dependence 
of the final radius on a is rather weak and non-monotonic because it competes with the 
dependence on the surface density whose precise value is an intricate function of a. On 

^ Indeed, we have the energy balance: oc Si/tf2| where To is the effective temperature, E the 
gas surface density, Ut the turbulent viscosity and the Keplerian angular velocity. Besides, 
the mid-plane temperature Tc obeys the relation T* ~ rT* where r is the opacity. On the other 
hand, the opacity r is given by r ~ kT, where k is the Rosseland opacity, while the turbulent 
viscosity scales as ft ~ aTc/flk, hence we obtain: oc naT,^ . 
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the other hand, the migration time usually increases going from a = 10~* up to a 0.1. 
Indeed, as seen from Fig. ^ a smaller a can lead to a larger mid-plane solid density. 
This is partly due to the fact that the dust disk height is smaller because the turbulence 
measured by a is weaker. 

The results displayed in Fig. |21 show that for a ^ 0.01 a Jupiter-like planet can 
migrate to a very small distance from the parent star, 0.03AU < r < O.lAU, provided 
the disk mass is sufficiently large Md ^ 10^"^ A/0. Only in the cases Md ^ IO^^A/q, or 
Md ^ lO~^Af0 with a ^ 0.1, the interaction with the planetesimal disk is too weak to 
yield a significant migration. Then, our results show that a final radius of 0.03AU < r < 
O.lAU is a natural result, unless the planetesimal disk has been cleared off too early on 
(e.g., by gravitational scattering). 

In summary, the present model predicts that, unless the disk mass is very small 
Md ^ IO^'^A/q, planets tend to move close to the parent star and to pile up to distances 
of the order of 0.03 — 0.04 AU (see Kuchner & Lecar 2002). However, with some degree 
of fine-tuning it is also possible to find a planet at intermediate distances between its 
formation site and such small radii. 

5. Conclusions 

In the current paper, we further improved the model for the migration of planets intro- 
duced in DPI and extended to time-dependent planetesimal accretion disks in DP2. After 
releasing the assumption of DP2 that the surface density of planetesimals is proportional 
to that of gas, a simplified model developed by Stepinski & Valageas (1996, 1997), that 
is able to simultaneously follow the evolution of gas and high-temperature silicates for 
up to lO^yr, is used. Then we coupled this disk model to the migration model introduced 
in DPI in order to obtain the migration rate of the planets in the planetesimal disk and 
to study how the migration rate depends on the disk mass, on its time evolution and on 
the dimensionless viscosity parameter a. 

We found that in the case of disks having a total mass of Md > IO^^Mq planets can 
migrate inward over a large distance, while if Af < lO~'^A/0 the planets remain almost 
at their initial position. On the other hand, for Md ^ lO^'^A'/0 a significant migration 
requires a ^ 10^^. 

If the migration is efficient the planet usually ends up at a small radius in the range 
0.03 — 0.1 AU which is simply set by the evaporation radius of the gaseous disk which 
gave rise at earlier times to the radial distribution of the planetesimal swarm. Therefore, 
our model provides a natural explanation for the small observed radii of extra-solar giant 
planets. In order to inhibit this process (so that Jupiter-like planets like our own remain 
at larger distances ^ 5 AU) the planetesimal disk must be cleared off over a time-scale of 
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Fig. 2. (a) The evolution of the semi-major axis a{t) of a Jupiter-mass planet, M = IMj, 
in a planctesimal disk for a = 0.1 and several values of M^, 0.1 (solid line), 0.01 (dotted 
line), 0.001 (short-dashed line) and O.OOOIA'/q (long-dashed line), (b) Same as Fig. (la) 
but with a — 0.01. (c) Same as Fig. (la) but with a = 0.001. (d) Same as Fig. (la) but 
with a = 0.0001. 

the order of or smaller than 10^ yr (depending on the properties of the disk) or the disk 
mass must be rather small (i.e. smaller than IMj) which could suggest an alternative 
formation scenario for such a giant planet (i.e. not related with the disk itself). 
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